Question 1: How do environmental conditions change by a) canopy treatment and b) canopy cover?

Is canopy treatment or canopy cover a stronger predictor of the environmental conditions below the shrub canopy?

Notes: * response variable: The variable you want to test. (e.g.,. PAR, TDR, nncover) * explanatory variables: The variables that can explain the outcome. (e.g., treatment, mafacover, site) * data frame: Include all the data you’ll be working with. (e.g., combine par and mafa_cover_all) * formula: a formula is a two-sided linear formula object describing both the fixed-effects and random effects part of the model, with the response on the left of a ~ operator and the terms, separated by + operators, on the right. * random effects: Random-effects terms are distinguished by vertical bars (|) separating expressions for design matrices from grouping factors. Two vertical bars (||) can be used to specify multiple uncorrelated random effects for the same grouping variable. + Because of the way it is implemented, the ||-syntax works only for design matrices containing numeric (continuous) predictors; to fit models with independent categorical effects, see dummy or the lmer_alt function from the afex package.)

Load PAR, TDR, and non-native percent cover data.

Save as a csv in the “Processed” data folder.

PAR linear mixed-effects model

  1. lmer
  2. testing the residuals for normality
lme.parmod <- lmer(normalized_par ~ treatment + mafacover + site + (1|date) +(1|plot_number),
                   data = parmod)

MOD <- lme.parmod

# Calculate model's residual and fitted values
F1 <- fitted(MOD)
E1 <- residuals(MOD, "pearson")

# See residuals as a histogram, look for normal distribution
hist(E1)

# Graph residuals against a normal distribution, look for all points to be on 1:1 line
#need to call next two lines together
qqnorm(E1, xlab = "Theoretical Quantiles", ylab = "Pearson Residuals")
qqline(E1)

# Graph residuals vs fitted values to inspect homogeneity of variance, look for a random scattering of plots
# run together 2 lines
plot(x=F1, y=E1, xlab="Fitted Model Values", ylab="Pearson Residuals", main = "Variance")
abline(h = 0) # looking for random scatter

# Check for crazy outliers, look for lines MUCH taller than the rest of data
plot(cooks.distance(MOD))

# Plot a box plot for each factor, look for equal variance amount levels along the zero line
# Box plots need to be assessed for every fixed effect (categorical factor) for your model. 
# Random effects need the qqline just like normality. 
# So, only site and treatment get box plots. 
# MAFA cover is a covariate (continuous variable).
boxplot(E1 ~ treatment, data = parmod)
abline(h = 0)

boxplot(E1 ~ site, data = parmod)
abline(h = 0)

# Plot each random effect in the model - things in the (1|_), look for dots along 1:1 line
qqnorm(ranef(MOD)$date[,1], main = "Plot Residuals")
qqline(ranef(MOD)$date[,1]) # unique dates

qqnorm(ranef(MOD)$plot_number[,1], main = "Plot Residuals")
qqline(ranef(MOD)$plot_number[,1]) # unique plots

3. statistics

# print 
print(lme.parmod)
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: normalized_par ~ treatment + mafacover + site + (1 | date) +  
##     (1 | plot_number)
##    Data: parmod
## REML criterion at convergence: -247.8414
## Random effects:
##  Groups      Name        Std.Dev.
##  plot_number (Intercept) 0.09588 
##  date        (Intercept) 0.08129 
##  Residual                0.13089 
## Number of obs: 288, groups:  plot_number, 36; date, 8
## Fixed Effects:
##   (Intercept)  treatmenthalf    treatmentno      mafacover           site  
##      0.413527      -0.004588      -0.011217      -0.003168       0.035283
summary(lme.parmod) # p-values
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: normalized_par ~ treatment + mafacover + site + (1 | date) +  
##     (1 | plot_number)
##    Data: parmod
## 
## REML criterion at convergence: -247.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9319 -0.7023 -0.0626  0.6010  3.1787 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  plot_number (Intercept) 0.009192 0.09588 
##  date        (Intercept) 0.006608 0.08129 
##  Residual                0.017131 0.13089 
## Number of obs: 288, groups:  plot_number, 36; date, 8
## 
## Fixed effects:
##                 Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)    4.135e-01  8.065e-02  4.954e+01   5.128 4.87e-06 ***
## treatmenthalf -4.588e-03  4.610e-02  3.458e+01  -0.100  0.92129    
## treatmentno   -1.122e-02  4.892e-02  3.856e+01  -0.229  0.81987    
## mafacover     -3.168e-03  9.918e-04  1.116e+02  -3.194  0.00183 ** 
## site           3.528e-02  4.077e-02  3.995e+01   0.866  0.39193    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trtmnth trtmntn mafcvr
## treatmnthlf -0.074                       
## treatmentno  0.009  0.572                
## mafacover   -0.541 -0.333  -0.459        
## site        -0.841 -0.164  -0.226   0.492
anova(lme.parmod) # fixed level effects, pvalues, etc.
## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq  Mean Sq NumDF   DenDF F value   Pr(>F)   
## treatment 0.000926 0.000463     2  34.385  0.0270 0.973361   
## mafacover 0.174745 0.174745     1 111.567 10.2005 0.001825 **
## site      0.012833 0.012833     1  39.948  0.7491 0.391927   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  1. figures
# scatterplot: par x mafacover
ggplot(data = parmod, aes(x = mafacover, y = normalized_par)) + 
  geom_point() +
  geom_smooth(method="lm")+
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

# scatterplot: par x site
ggplot(data = parmod, aes(x = site, y = normalized_par)) + 
  geom_boxplot(aes(fill = site)) +
  theme_bw()
## Warning: Continuous x aesthetic
## ℹ did you forget `aes(group = ...)`?
## Warning: The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

# scatterplot: mafacover x par x treatment
ggplot(data = parmod, aes(x = mafacover, y = normalized_par)) + 
  geom_point(aes(color = treatment)) +
  theme_bw()

# boxplot: mafacover x site
ggplot(data = parmod, aes(x = site, y = normalized_par)) +
  geom_boxplot(aes(fill = site)) +
  theme_bw() + 
  theme(legend.position = "none")
## Warning: Continuous x aesthetic
## ℹ did you forget `aes(group = ...)`?
## The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

parmod_treatment <- aov(normalized_par ~ treatment, data = parmod) 
TukeyHSD(parmod_treatment)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = normalized_par ~ treatment, data = parmod)
## 
## $treatment
##                  diff         lwr          upr     p adj
## half-full -0.05365265 -0.11695323  0.009647926 0.1147974
## no-full   -0.08296590 -0.14626648 -0.019665321 0.0062500
## no-half   -0.02931325 -0.09261383  0.033987332 0.5204084
plot(TukeyHSD(parmod_treatment))

TDR linear mixed-effects model

  1. lmer
  2. testing the residuals for normality
lme.tdrmod <- lmer(soil_moisture ~ treatment + mafacover + site + (1|date) +(1|plot_number),
                   data = tdrmod)

MOD <- lme.tdrmod

# Calculate model's residual and fitted values
F1 <- fitted(MOD)
E1 <- residuals(MOD, "pearson")

# See residuals as a histogram, look for normal distribution
hist(E1)

# Graph residuals against a normal distribution, look for all points to be on 1:1 line
#need to call next two lines together
qqnorm(E1, xlab = "Theoretical Quantiles", ylab = "Pearson Residuals")
qqline(E1)

# Graph residuals vs fitted values to inspect homogeneity of variance, look for a random scattering of plots
# run together 2 lines
plot(x=F1, y=E1, xlab="Fitted Model Values", ylab="Pearson Residuals", main = "Variance")
abline(h = 0) # looking for random scatter

# Check for crazy outliers, look for lines MUCH taller than the rest of data
plot(cooks.distance(MOD))

# Plot a box plot for each factor, look for equal variance amount levels along the zero line
# Box plots need to be assessed for every fixed effect (categorical factor) for your model. 
# Random effects need the qqline just like normality. 
# So, only site and treatment get box plots. 
# MAFA cover is a covariate (continuous variable).
boxplot(E1 ~ treatment, data = tdrmod)
abline(h = 0)

boxplot(E1 ~ site, data = tdrmod)
abline(h = 0)

# Plot each random effect in the model, look for dots along 1:1 line
qqnorm(ranef(MOD)$date[,1], main = "Plot Residuals")
qqline(ranef(MOD)$date[,1])

qqnorm(ranef(MOD)$plot_number[,1], main = "Plot Residuals")
qqline(ranef(MOD)$plot_number[,1])

3. statistics

# print 
print(lme.tdrmod)
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: soil_moisture ~ treatment + mafacover + site + (1 | date) + (1 |  
##     plot_number)
##    Data: tdrmod
## REML criterion at convergence: 522.5547
## Random effects:
##  Groups      Name        Std.Dev.
##  plot_number (Intercept) 1.372   
##  date        (Intercept) 7.207   
##  Residual                2.249   
## Number of obs: 108, groups:  plot_number, 36; date, 4
## Fixed Effects:
##   (Intercept)  treatmenthalf  treatmentnone      mafacover           site  
##      21.81315       -0.86322       -3.15994        0.01808       -4.06263
summary(lme.tdrmod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: soil_moisture ~ treatment + mafacover + site + (1 | date) + (1 |  
##     plot_number)
##    Data: tdrmod
## 
## REML criterion at convergence: 522.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.07912 -0.53157  0.00198  0.51062  2.23084 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  plot_number (Intercept)  1.884   1.372   
##  date        (Intercept) 51.934   7.207   
##  Residual                 5.060   2.249   
## Number of obs: 108, groups:  plot_number, 36; date, 4
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   21.81315    3.90031  4.07414   5.593  0.00475 ** 
## treatmenthalf -0.86322    0.82965 35.78406  -1.040  0.30511    
## treatmentnone -3.15994    0.92207 41.74189  -3.427  0.00138 ** 
## mafacover      0.01808    0.02545 79.56419   0.710  0.47961    
## site          -4.06263    0.83253 55.27825  -4.880  9.4e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trtmnth trtmntn mafcvr
## treatmnthlf -0.006                       
## treatmentnn  0.045  0.591                
## mafacover   -0.233 -0.368  -0.548        
## site        -0.352 -0.199  -0.297   0.542
anova(lme.tdrmod) # A significant p-value indicates that 
## Type III Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## treatment  66.942  33.471     2 36.100  6.6149  0.003568 ** 
## mafacover   2.553   2.553     1 79.564  0.5045  0.479615    
## site      120.493 120.493     1 55.278 23.8131 9.403e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# one or more significant differences exist between group means.

tdrmod_site <- aov(soil_moisture ~ site*treatment, data = tdrmod) 
anova(tdrmod_site)
## Analysis of Variance Table
## 
## Response: soil_moisture
##                 Df Sum Sq Mean Sq F value  Pr(>F)   
## site             1  344.6  344.57  7.7568 0.00638 **
## treatment        2  154.9   77.45  1.7434 0.18009   
## site:treatment   2    4.6    2.29  0.0516 0.94971   
## Residuals      102 4531.0   44.42                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(tdrmod_site)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: site
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: site,
## treatment
## Warning in TukeyHSD.aov(tdrmod_site): 'which' specified some non-factors which
## will be dropped
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = soil_moisture ~ site * treatment, data = tdrmod)
## 
## $treatment
##                 diff       lwr      upr     p adj
## half-full -0.6462963 -4.382655 3.090062 0.9110153
## none-full -2.8011574 -6.537516 0.935201 0.1803794
## none-half -2.1548611 -5.891219 1.581497 0.3595155
plot(TukeyHSD(tdrmod_site))
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: site
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: site,
## treatment
## Warning in TukeyHSD.aov(tdrmod_site): 'which' specified some non-factors which
## will be dropped

4. figures

# boxplot: TDR by treatment by site
ggplot(data = tdrmod, aes(x = treatment, y = soil_moisture)) +
  geom_boxplot(aes(fill = site)) +
  theme_bw() 
## Warning: The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

# boxplot: TDR by site by treatment
ggplot(data = tdrmod, aes(x = site, y = soil_moisture)) +
  geom_boxplot(aes(fill = treatment)) +
  scale_fill_manual(values = c('#F6EECF', '#B09175', '#291611')) +
  labs(title = "Soil Moisture by Site and Treatment", x = "Site", y = "PAR") +
  theme_bw() 

# scatterplot: tdr x mafacover
ggplot(data = tdrmod, aes(x = mafacover, y = soil_moisture)) + 
  geom_point() +
  geom_smooth(method="lm")+
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

Non-native cover linear mixed-effects model

  1. lmer
  2. testing the residuals for normality
#need to convert df to a tibble
#as_tibble(nnmod)
#print(nn)
#names(nnmod)
#print(nnmod)
#str(nnmod)
lme.nnmod <- lmer(nncover ~ treatment + mafacover + site + (1|date) + (1|plot_number),
                   data = nnmod)

MOD <- lme.nnmod

# Calculate model's residual and fitted values
F1 <- fitted(MOD)
E1 <- residuals(MOD, "pearson")

# See residuals as a histogram, look for normal distribution
hist(E1)

# Graph residuals against a normal distribution, look for all points to be on 1:1 line
#need to call next two lines together
qqnorm(E1, xlab = "Theoretical Quantiles", ylab = "Pearson Residuals")
qqline(E1)

# Graph residuals vs fitted values to inspect homogeneity of variance, look for a random scattering of plots
# run together 2 lines
plot(x=F1, y=E1, xlab="Fitted Model Values", ylab="Pearson Residuals", main = "Variance")
abline(h = 0) # looking for random scatter

# Check for crazy outliers, look for lines MUCH taller than the rest of data
plot(cooks.distance(MOD))

# Plot a box plot for each factor, look for equal variance amount levels along the zero line
# Box plots need to be assessed for every fixed effect (categorical factor) for your model. 
# Random effects need the qqline just like normality. 
# So, only site and treatment get box plots. 
# MAFA cover is a covariate (continuous variable).
boxplot(E1 ~ treatment, data = nnmod)
abline(h = 0)

boxplot(E1 ~ site, data = nnmod)
abline(h = 0)

# Plot each random effect in the model, look for dots along 1:1 line
qqnorm(ranef(MOD)$date[,1], main = "Plot Residuals")
qqline(ranef(MOD)$date[,1])

qqnorm(ranef(MOD)$plot_number[,1], main = "Plot Residuals")
qqline(ranef(MOD)$plot_number[,1])

3. statistics

  1. figures
# boxplot: nncover x treatment x site
ggplot(data = nnmod, aes(x = treatment, y = nncover)) + 
  geom_boxplot(aes(fill = site)) +
  geom_smooth(method="lm") +
  labs(title = "Non-native cover by Treatment and Site", x = "Treatment", y = "Non-native cover (%)") +
  scale_fill_manual(values = c( "#DCC27A", "#F19B34")) +
  # scale_fill_manual(values = c('#D3E3CA', '#92A587')) +
  theme_bw()
## Warning: The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## `geom_smooth()` using formula = 'y ~ x'

# boxplot: mafacover x site x treatment
ggplot(data = nnmod, aes(x = site, y = nncover)) + 
  geom_boxplot(aes(fill = treatment)) +
 # scale_fill_manual(values = c('#F19B34', '#9F7E75', '#92A587')) +
  scale_fill_manual(values = c('#F6EECF', '#B09175', '#291611')) +
  labs(title = "Non-native cover by Site and Treatment", x = "Site", y = "Non-native cover (%)") +
  geom_smooth(method="lm") +
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

# scatterplot: nncover x mafacover
ggplot(data = nnmod, aes(x = mafacover, y = nncover)) + 
  geom_point() +
  geom_smooth(method="lm")+
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

Take aways

#Question 2: How does seedling survival change by PAR, TDR, and mafacover? Is canopy treatment or canopy cover a stronger predictor of the environmental conditions below the shrub canopy?

Glossary

df exploration * names(df) will give you the column names in the df

statistics * lmer: Fit Linear Mixed-Effects Models + Fit a linear mixed-effects model (LMM) to data, via REML or maximum likelihood.

end of code